//
// Classe de convolucao convolve.cpp
//

#ifndef __CONVOLVE__
#define __CONVOLVE__

#include <iostream>
#include <gsl/gsl_interp.h>
#include <gsl/gsl_integration.h>

class GSL_Integration
{
	gsl_integration_workspace *w;
	gsl_integration_cquad_workspace *_wcquad;
	gsl_function *F;
	double limits[2];
	double res;
	double err;
	bool _integrated;
	size_t _size;
	
public:
	GSL_Integration(gsl_function &func,int size = 1000)
	{
	w = gsl_integration_workspace_alloc(size);
	_wcquad = gsl_integration_cquad_workspace_alloc(size);
	F=&func;
	_size=size;
	_integrated=false;
	};
	~GSL_Integration(void)
	{
#ifdef DEBUG
	cout << "# - GSL_Integration destructor called...\n";
#endif
	gsl_integration_workspace_free (w);
	gsl_integration_cquad_workspace_free (_wcquad);
	
	};
	double integrate_qag(void)
	{
	gsl_integration_qags (F, limits[0], limits[1], 1e2, 1e-1, _size,w, &res, &err);
		//cout << "aqui\n";
	_integrated=true;
	return res;
	};
	double integrate_qng(void)
	{
	gsl_integration_qng (F, limits[0], limits[1],1e-1, 1e-1, &res, &err, &_size);
	_integrated=true;
	return res;
	
	};
	double integrate_cquad(void)
	{
	gsl_integration_cquad (F, limits[0], limits[1],1e-1, 1e-1, _wcquad, &res,&err,&_size);
	_integrated=true;
	return res;
	
	}
	void inline clear(void)
	{
	gsl_integration_workspace_free (w);
	w = gsl_integration_workspace_alloc(_size);
	gsl_integration_cquad_workspace_free (_wcquad);
	_wcquad = gsl_integration_cquad_workspace_alloc(_size);
	
	};
	void inline clean(void)
	{
	gsl_integration_workspace_free (w);
	gsl_integration_cquad_workspace_free (_wcquad);
	};
	void xmin(double val){limits[0]=val;_integrated=false;};
	void xmax(double val){limits[1]=val;_integrated=false;};
	double xmin(void){return limits[0];};
	double xmax(void){return limits[1];};
	double inline integrationErr(void){return err;};
	
};



//using namespace std;

template <class T> 
class GSL_Interp{
	double *x1_axis;
	double *y1_axis;
	double _integr1;
	int _sz1;
	gsl_interp *f1_interp;
	gsl_interp_accel *acc1;
	double _a;
	double _b;
	public:

GSL_Interp(T &f1,const int size1)
	{
		//T = f1;
		f1_interp = gsl_interp_alloc(gsl_interp_linear,size1);
		x1_axis = new double [size1];
		y1_axis = new double [size1];
		_sz1=size1;
		for (int i = 0; i < size1; i++){
			x1_axis[i]=f1[0][i];
			y1_axis[i]=f1[1][i];
		}
		acc1 = gsl_interp_accel_alloc ();
		gsl_interp_init(f1_interp,x1_axis,y1_axis,size1);
		_integr1 = gsl_interp_eval_integ(f1_interp,x1_axis,y1_axis,f1[0][0],f1[0][_sz1-1],acc1);
		_a = f1[0][0];
		_b = f1[0][_sz1-1];
	};
	~GSL_Interp()
	{
		if(f1_interp) gsl_interp_free(f1_interp);
		if(acc1) gsl_interp_accel_free (acc1);
		if(x1_axis) delete [] x1_axis;
		if(y1_axis) delete [] y1_axis;

	};
	double inline f(const double px)
	{
		return gsl_interp_eval(f1_interp,x1_axis,y1_axis,px,acc1);
	};
	double inline integrate(void){return _integr1;};
	double inline xmin(void){return _a;};
	double inline xmax(void){return _b;};
	void inline xmin(double val)
	{
		_a=val;
		_integr1 = gsl_interp_eval_integ(f1_interp,x1_axis,y1_axis,_a,_b,acc1);
	};
	void inline xmax(double val)
	{
		_b=val;
		_integr1 = gsl_interp_eval_integ(f1_interp,x1_axis,y1_axis,_a,_b,acc1);
	};
};

class Convolve{
		double *x1_axis;
		double *y1_axis;
		double *x2_axis;
		double *y2_axis;
		double _integr1;
		double _integr2;
		int _sz1;
		int _sz2;
		gsl_interp *f1_interp;
		gsl_interp *f2_interp;
		gsl_interp_accel *acc1;
		gsl_interp_accel *acc2;
		double _px;
		double _a;
		double _b;
	public:
		Convolve(double **f1,double **f2, const int size1, const int size2)
		{
		//Achar funcao que descreve x e y -> interpolacao linear
			//_px=0.0;
			f1_interp = gsl_interp_alloc(gsl_interp_linear,size1);
			f2_interp = gsl_interp_alloc(gsl_interp_linear,size2);
			x1_axis = new double [size1];
			y1_axis = new double [size1];
			x2_axis = new double [size2];
			y2_axis = new double [size2];
			_sz1 = size1;
			_sz2 = size2;
			for (int i = 0; i < size1; i++){
				x1_axis[i]=f1[0][i];
				y1_axis[i]=f1[1][i];
			}
			for (int i = 0; i < size2; i++){
				x2_axis[i]=f2[0][i];
				y2_axis[i]=f2[1][i];
			}
			acc1 = gsl_interp_accel_alloc ();
			acc2 = gsl_interp_accel_alloc ();
			gsl_interp_init(f1_interp,f1[0],f1[1],size1);
			gsl_interp_init(f2_interp,f2[0],f2[1],size2);
			_integr1 = gsl_interp_eval_integ(f1_interp,f1[0],f1[1],f1[0][0],f1[0][_sz1-1],acc1);
			_integr2 = gsl_interp_eval_integ(f2_interp,f2[0],f2[1],f2[0][0],f2[0][_sz1-1],acc2);
			_px=0.0;
			//cout << "# "<< _integr1 << "  " << _integr2 << endl;
		};
		Convolve(vector<vector<double> > &f1,vector<vector<double> > &f2, const int size1, const int size2)
		{
		//Achar funcao que descreve x e y -> interpolacao linear
			//_px=0.0;
			f1_interp = gsl_interp_alloc(gsl_interp_linear,size1);
			f2_interp = gsl_interp_alloc(gsl_interp_linear,size2);
			x1_axis = new double [size1];
			y1_axis = new double [size1];
			x2_axis = new double [size2];
			y2_axis = new double [size2];
			_sz1 = size1;
			_sz2 = size2;
			for (int i = 0; i < size1; i++){
				x1_axis[i]=f1[0][i];
				y1_axis[i]=f1[1][i];
			}
			for (int i = 0; i < size2; i++){
				x2_axis[i]=f2[0][i];
				y2_axis[i]=f2[1][i];
			}
			acc1 = gsl_interp_accel_alloc ();
			acc2 = gsl_interp_accel_alloc ();
			gsl_interp_init(f1_interp,x1_axis,y1_axis,size1);
			gsl_interp_init(f2_interp,x2_axis,y2_axis,size2);
			_integr1 = gsl_interp_eval_integ(f1_interp,x1_axis,y1_axis,x1_axis[0],x1_axis[_sz1-1],acc1);
			_integr2 = gsl_interp_eval_integ(f2_interp,x2_axis,x2_axis,x2_axis[0],x2_axis[_sz1-1],acc2);
			_px=0.0;
			//cout << _integr1 << "  " << _integr2 << endl;
		};
	
		~Convolve(void)
		{/*
			if(f1_interp) gsl_interp_free(f1_interp);
			if(f2_interp) gsl_interp_free(f2_interp);
			if(acc1) gsl_interp_accel_free (acc1);
			if(acc2) gsl_interp_accel_free (acc2);
			if(x1_axis) delete [] x1_axis;
			if(y1_axis) delete [] y1_axis;
			if(x2_axis) delete [] x2_axis;
			if(y2_axis) delete [] y2_axis;
			*/
			//gsl_interp_free(f1_interp);
			//gsl_interp_free(f2_interp);
			//gsl_interp_accel_free (acc1);
			//gsl_interp_accel_free (acc2);
			//delete [] x1_axis;
			//delete [] y1_axis;
			//delete [] x2_axis;
			//delete [] y2_axis;
			
			//gsl_interp_free(f1_interp);
			//gsl_interp_free(f2_interp);
			//delete [] x1_axis;
			//delete [] y1_axis;
			//delete [] x2_axis;
			//delete [] y2_axis; 
		};
	/*	  */
	Convolve(double **f1,double (f2)(double x), const int size1)
	{
		//Achar funcao que descreve x e y -> interpolacao linear
		//_px=0.0;
		f1_interp = gsl_interp_alloc(gsl_interp_linear,size1);
		x1_axis = new double [size1];
		y1_axis = new double [size1];
		_sz1 = size1;
		for (int i = 0; i < size1; i++){
			x1_axis[i]=f1[0][i];
			y1_axis[i]=f1[1][i];
		}
		acc1 = gsl_interp_accel_alloc ();
		gsl_interp_init(f1_interp,f1[0],f1[1],size1);
		_integr1 = gsl_interp_eval_integ(f1_interp,f1[0],f1[1],f1[0][0],f1[0][_sz1-1],acc1);
		_px=0.0;
	};
	
		void inline limits(double a, double b)
		{
			_a=a;
			_b=b;
		};
		double inline integrate1(void){return _integr1;};			
		double inline integrate2(void){return _integr2;};
		double inline xmin(void){return _a;};
		double inline xmax(void){return _b;};
		void inline xmin(double a){_a=a;};
		void inline xmax(double b){_b=b;};
		void inline px(double val){_px=val;};
		double inline px(void){return _px;};
		void inline clean(void){
			
			gsl_interp_free(f1_interp);
			gsl_interp_free(f2_interp);
			gsl_interp_accel_free (acc1);
			gsl_interp_accel_free (acc2);
			delete [] x1_axis;
			delete [] y1_axis;
			delete [] x2_axis;
			delete [] y2_axis;
			 
		};
		double inline cnv_fnc(double x, double t)
		{
			double px = t-x; // *(double *) params;
			//cout << "aqui 0 " << x << endl; 
			double a = 0.0; 
			if((x1_axis[0] < x) && (x < x1_axis[_sz1-1]))
				a = gsl_interp_eval(f1_interp,x1_axis,y1_axis,x,acc1);
			//cout << "aqui 1\n";
			double b = 0.0;
			if((xmin() < px) && (px < xmax())){
				//cout << "aqui\n";
				b = gsl_interp_eval(f2_interp,x2_axis,y2_axis,px,acc2);
			}
			//cout << "aqui 2\n";
			return a*b; // *_integr2);
		};

//		double inline val(const double px)
//		{
//			Integrate(&cnv_fnc,px) res;
//			return res.px;
//		};

};

double cnv_fnc(double x, void *par)
{
	Convolve ff = *( Convolve *) par;
	double val = ff.cnv_fnc(x,ff.px());
		//cout << " -> " << x << " " << ff.px() << endl;
		//cout  << x << " " << val << endl;
	return val;
}

#endif